      subroutine control(spiu, smeno, cpiu, cmeno) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use corgan_com_M 
      use cindex_com_M 
      use numpar_com_M 
      use cprplt_com_M, ONLY:  
      use cophys_com_M 
      use ctemp_com_M, ONLY: fmf, fef, fvf, fbxf, fbyf, fbzf, twx 
      use blcom_com_M, ONLY: npcel_control, lsplit, lcoal, mtresh, vtresh, &
         number, x, y, z, iphead, px, py, pz, pxi, peta, pzta, up, vp, wp, &
         link, ico, qpar 
      use objects_com_M, ONLY:  
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
!
!     Library of routines to control the number
!     of particles in 3d kinetic plasma
!     simulations -Derived from the 2d version for Celeste2d
!
!     by Gianni Lapenta
!
 
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      real(double)  :: spiu 
      real(double)  :: smeno 
      real(double)  :: cpiu 
      real(double)  :: cmeno 
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv 
      integer , dimension(4) :: lioinput, lioprint 
      integer :: ncount, np, is, icoal, ncounta 
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2 
      real(double) :: ixi4, ieta4, ixi5, ieta5 
      real(double), dimension(8) :: wght 
!      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
!         ucm, vcm, wcm, cm 
!      real(double), dimension(3,20) :: wsin, wcos 
!      real(double), dimension(3,12,20) :: tsi, rot 
      real(double) :: cmenold 
      logical :: expnd, nomore, testo 
!      character :: name*80 
!-----------------------------------------------
      testo = .true. 
!
      kbp1 = kbp2 - 1 
      jbp1 = jbp2 - 1 
      ibp1 = ibp2 - 1 
!
      if (testo) then 
         ncount = 0 
         np = iphead(1) 
   10    continue 
         if (np /= 0) then 
            ncount = ncount + 1 
            np = link(np) 
            go to 10 
         endif 
         write (6, *) 'avail. particles', ncount 
      endif 
 
      do is = 1, nsp 
         do k = 2, kbp1 
            do j = 2, jbp1 
               do i = 2, ibp1 
                  ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
                  if (number(ijk,is)<npcel_control(is) - 5 .and. lsplit) call &
                     splitter (ijk, i, j, k, is, spiu, smeno, iphead, link, &
                     wght, up, vp, wp, pxi, peta, pzta, px, py, pz, qpar, ico, &
                     iwid, jwid, kwid, x, y, z, npcel_control) 
                  ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
                  if (.not.(number(ijk,is)>npcel_control(is)+5. .and. lcoal)) &
                     cycle  
!     write(6,*)'coalesce called',cpiu,cmeno,number(ij,is),ij,is
                  do icoal = 1, 1!int(number(ijk,is)-npcel_control(is)) + 1 
                     cmenold = cmeno 
                     if (testo) then 
                        ncounta = 0 
                        np = iphead(ijk) 
    1                   continue 
                        if (np /= 0) then 
                           ncounta = ncounta + 1 
                           np = link(np) 
                           go to 1 
                        endif 
                     endif 
                     call coalesce (ijk, i, j, k, is, cpiu, cmeno, mtresh, &
                        vtresh, iphead, link, wght, up, vp, wp, pxi, peta, pzta&
                        , px, py, pz, qpar, ico, iwid, jwid, kwid, x, y, z) 
                     if (testo) then 
                        ncount = 0 
                        np = iphead(ijk) 
   11                   continue 
                        if (np /= 0) then 
                           ncount = ncount + 1 
                           np = link(np) 
                           go to 11 
                        endif 
                        if (ncounta < ncount) write (*, *) 'coalesce error  ', &
                           ncount, ncounta, ijk, is, cmeno - cmenold 
                     endif 
                     if (cmeno <= cmenold) cycle  
                     exit  
                  end do 
               end do 
            end do 
         end do 
      end do 
      return  
      end subroutine control 


 
      subroutine coalesce(ijk, i, j, k, is, suc, fail, ntresh, treshco, iphead&
         , link, wght, up, vp, wp, xp, yp, zp, xpf, ypf, zpf, qpar, ico, iwid, &
         jwid, kwid, x, y, z) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
!
!     A routine to coalesce 2 particles in 3d kinetic
!     problems
!
      USE vast_kind_param, ONLY:  double 
      use controllo_M
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ijk 
      integer  :: i 
      integer  :: j 
      integer  :: k 
      integer , intent(in) :: is 
      integer , intent(in) :: ntresh 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      real(double) , intent(inout) :: suc 
      real(double) , intent(inout) :: fail 
      real(double) , intent(in) :: treshco 
      integer , intent(inout) :: iphead(*) 
      integer  :: link(0:*) 
      integer , intent(inout) :: ico(0:*) 
      real(double)  :: wght(*) 
      real(double) , intent(inout) :: up(0:*) 
      real(double) , intent(inout) :: vp(0:*) 
      real(double) , intent(inout) :: wp(0:*) 
      real(double) , intent(inout) :: xp(0:*) 
      real(double) , intent(inout) :: yp(0:*) 
      real(double) , intent(inout) :: zp(0:*) 
      real(double) , intent(out) :: xpf(0:*) 
      real(double) , intent(out) :: ypf(0:*) 
      real(double) , intent(out) :: zpf(0:*) 
      real(double) , intent(inout) :: qpar(0:*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      real(double) :: e1, e2, xl, xr, xm, yb, yt, ym, ze, zf, zm, xpp, ypp, zpp, dmmin&
         , qtotl, ume, u2me, dist, vme, v2me, wme, w2me, xi, eta, zta 
      logical :: flag 
!-----------------------------------------------
!   E x t e r n a l   F u n c t i o n s
!-----------------------------------------------
      integer , external :: cvmgp 
!-----------------------------------------------
!
!
!
      fail = fail + 1 
      e1 = 1.E10 
      e2 = 1.E10 
      np1 = 0 
      ipdone = 0 
      np2 = 0 
!
      np = iphead(ijk) 
      if (iphead(ijk) /= 0) then 
 
!
!     dicotomic search
!
 
 
         call listmkr (np, nploc, link, lvmax, ntmp) 
 
         flag = .TRUE. 
         xl = 0. 
         xr = 1. 
         xm = 0.5 
         yb = 0. 
         yt = 1. 
         ym = 0.5 
         ze = 0. 
         zf = 1. 
         zm = 0.5 
         ntimes = 0 
 
   11    continue 
         ntimes = ntimes + 1 
         nvlist = 0 
!DIR$ IVDEP
         do lv = 1, lvmax 
            np = nploc(lv) 
            if (ico(np) /= is) cycle  
            xpp = xp(np) - int(xp(np)) 
            ypp = yp(np) - int(yp(np)) 
            zpp = zp(np) - int(zp(np)) 
            ix = cvmgp(0,1,xm - xpp) 
            iy = cvmgp(0,1,ym - ypp) 
            iz = cvmgp(0,1,zm - zpp) 
            itmp(lv) = 1 + ix + iy*2 + iz*4 
         end do 
!
         do lv = 1, lvmax 
!
            np = nploc(lv) 
            if (ico(np) /= is) cycle  
            nvlist(itmp(lv)) = nvlist(itmp(lv)) + 1 
            lvlist(itmp(lv),nvlist(itmp(lv))) = np 
         end do 
!
         lvmax = 0 
         lismax = 0 
         do ii = 1, 8 
            if (nvlist(ii) < lvmax) cycle  
            lvmax = nvlist(ii) 
            lismax = ii 
         end do 
!DIR$ IVDEP
         nploc(:lvmax) = lvlist(lismax,:lvmax) 
         kmax = lismax/4 
         jmax = (lismax - 4*kmax)/2 
         imax = lismax - 4*kmax - 2*jmax 
         if (imax == 1) then 
            xr = xm 
         else 
            xl = xm 
         endif 
         if (jmax == 1) then 
            yt = ym 
         else 
            yb = ym 
         endif 
         if (kmax == 1) then 
            zf = zm 
         else 
            ze = zm 
         endif 
         xm = 0.5*(xl + xr) 
         ym = 0.5*(yb + yt) 
         zm = 0.5*(ze + zf) 
 
         if (lvmax>ntresh .and. ntimes<10) go to 11 
 
         if (lvmax == 2) then 
            np1 = nploc(1) 
            np2 = nploc(2) 
         else 
            dmmin = 1.E10 
            do lv1 = 1, lvmax - 1 
               n1 = nploc(lv1) 
               do lv2 = lv1 + 1, lvmax 
                  n2 = nploc(lv2) 
                  wdif(lv2) = sqrt((up(n1)-up(n2))**2+(vp(n1)-vp(n2))**2+(wp(n1&
                     )-wp(n2))**2) 
               end do 
!
               do lv2 = lv1 + 1, lvmax 
                  n2 = nploc(lv2) 
                  if (dmmin < wdif(lv2)) cycle  
                  dmmin = wdif(lv2) 
                  np1 = n1 
                  np2 = n2 
               end do 
            end do 
         endif 
 
 
         if (np1<=0 .or. np2<=0) then 
            write (66, *) 'coalesce failure in dicotomic search', np1, np2 
            return  
         endif 
 
 
 
         if (qpar(np1)*qpar(np2) < 0.) then 
            write (66, *) 'coalesce error: Particles of diff. charge' 
            return  
         endif 
         if (ico(np1)/=is .or. ico(np2)/=is) then 
            write (66, *) 'coalesce error: Particles of wrong species' 
            return  
         endif 
 
 
 
         qtotl = qpar(np1) + qpar(np2) 
         ume = (qpar(np1)*up(np1)+qpar(np2)*up(np2))/qtotl 
         u2me = sqrt((qpar(np1)*up(np1)**2+qpar(np2)*up(np2)**2)/qtotl) 
         dist = abs(abs(ume) - u2me)/(abs(ume) + u2me + 1E-10)*2 
         if (dist > treshco) then 
            write (66, *) 'coalesce fails the tests' 
            return  
         endif 
 
         vme = (qpar(np1)*vp(np1)+qpar(np2)*vp(np2))/qtotl 
         v2me = sqrt((qpar(np1)*vp(np1)**2+qpar(np2)*vp(np2)**2)/qtotl) 
         dist = abs(abs(vme) - v2me)/(abs(vme) + v2me + 1E-10)*2 
         if (dist > treshco) then 
            write (66, *) 'coalesce fails the tests' 
            return  
         endif 
 
         wme = (qpar(np1)*wp(np1)+qpar(np2)*wp(np2))/qtotl 
         w2me = sqrt((qpar(np1)*wp(np1)**2+qpar(np2)*wp(np2)**2)/qtotl) 
         dist = abs(abs(wme) - w2me)/(abs(wme) + w2me + 1E-10)*2 
         if (dist > treshco) then 
            write (66, *) 'coalesce fails the tests' 
            return  
         endif 
 
!
!
         np = iphead(ijk) 
   12    continue 
         if (np /= 0) then 
!
            if (np==np1 .or. np==np2) then 
!
               iphead(ijk) = link(np) 
               link(np) = iphead(1) 
               iphead(1) = np 
!
            else 
!
               iphead(ijk) = link(np) 
               link(np) = ipdone 
               ipdone = np 
!
            endif 
!
            np = iphead(ijk) 
!
            go to 12 
!
         endif 
!
!
         npn = iphead(1) 
         if (npn==0 .or. np1==0 .or. np2==0) then 
            write (*, *) 'coalesce major failure' 
            return  
         endif 
!
!     combine particles using simple linear combination
!
         xp(npn) = (xp(np1)*qpar(np1)+xp(np2)*qpar(np2))/qtotl 
         yp(npn) = (yp(np1)*qpar(np1)+yp(np2)*qpar(np2))/qtotl 
         zp(npn) = (zp(np1)*qpar(np1)+zp(np2)*qpar(np2))/qtotl 
!
         xi = xp(npn) - int(xp(npn)) 
         eta = yp(npn) - int(yp(npn)) 
         zta = zp(npn) - int(zp(npn)) 
!
         call weights (xi, eta, zta, wght) 
!
         ipjk = ijk + iwid 
         ipjpk = ijk + iwid + jwid 
         ijpk = ijk + jwid 
!
         ijkp = ijk + kwid 
         ipjkp = ijk + iwid + kwid 
         ijpkp = ijk + jwid + kwid 
         ipjpkp = ijk + iwid + jwid + kwid 
!
         xpf(npn) = wght(1)*x(ipjk) + (wght(2)*x(ipjpk)+(wght(3)*x(ijpk)+(wght(&
            4)*x(ijk)+(wght(5)*x(ipjkp)+(wght(6)*x(ipjpkp)+(wght(7)*x(ijpkp)+&
            wght(8)*x(ijkp))))))) 
!
         ypf(npn) = wght(1)*y(ipjk) + (wght(2)*y(ipjpk)+(wght(3)*y(ijpk)+(wght(&
            4)*y(ijk)+(wght(5)*y(ipjkp)+(wght(6)*y(ipjpkp)+(wght(7)*y(ijpkp)+&
            wght(8)*y(ijkp))))))) 
!
         zpf(npn) = wght(1)*z(ipjk) + (wght(2)*z(ipjpk)+(wght(3)*z(ijpk)+(wght(&
            4)*z(ijk)+(wght(5)*z(ipjkp)+(wght(6)*z(ipjpkp)+(wght(7)*z(ijpkp)+&
            wght(8)*z(ijkp))))))) 
!
!
!
!      up(npn)=sqrt(v2me)*vme/(1e-10+abs(vme))
         up(npn) = ume 
!      vp(npn)=sqrt(v2me)*vme/(1e-10+abs(vme))
         vp(npn) = vme 
!      wp(npn)=sqrt(v2me)*vme/(1e-10+abs(vme))
         wp(npn) = wme 
         qpar(npn) = qtotl 
         ico(npn) = is 
!
         iphead(1) = link(npn) 
         link(npn) = ipdone 
         ipdone = npn 
         suc = suc + 1 
         fail = fail - 1 
!
!
!
         iphead(ijk) = ipdone 
         ipdone = 0 
!
      endif 
!
      return  
      end subroutine coalesce 


 
      integer function isegno (x) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1   
     USE vast_kind_param, ONLY:  double             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      real(double) , intent(in) :: x 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
      if (x <= 0) then 
         isegno = 0
      else 
         isegno = 1
      endif 
      return  
      end function isegno 


 
 
 
      subroutine splitter(ijk, i, j, k, is, succ, fall, iphead, link, wght, up&
         , vp, wp, xp, yp, zp, xpf, ypf, zpf, qpar, ico, iwid, jwid, kwid, x, y&
         , z, npcel) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1    
      USE vast_kind_param, ONLY:  double          
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ijk 
      integer , intent(in) :: i 
      integer , intent(in) :: j 
      integer , intent(in) :: k 
      integer , intent(in) :: is 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      real(double) , intent(inout) :: succ 
      real(double) , intent(inout) :: fall 
      integer , intent(inout) :: iphead(*) 
      integer , intent(inout) :: link(0:*) 
      integer , intent(inout) :: ico(0:*) 
      integer , intent(in) :: npcel(*) 
      real(double)  :: wght(*) 
      real(double) , intent(inout) :: up(0:*) 
      real(double) , intent(inout) :: vp(0:*) 
      real(double) , intent(inout) :: wp(0:*) 
      real(double) , intent(inout) :: xp(0:*) 
      real(double) , intent(inout) :: yp(0:*) 
      real(double) , intent(inout) :: zp(0:*) 
      real(double) , intent(out) :: xpf(0:*) 
      real(double) , intent(out) :: ypf(0:*) 
      real(double) , intent(out) :: zpf(0:*) 
      real(double) , intent(inout) :: qpar(0:*) 
      real(double) , intent(in) :: x(*) 
      real(double) , intent(in) :: y(*) 
      real(double) , intent(in) :: z(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(6) :: npnew 
      integer :: nsplit, np, ip, ipjk, ipjpk, ijpk, ijkp, ipjkp, ijpkp, ipjpkp 
      real(double) , dimension(6) :: xpnew, ypnew, zpnew 
      real(double) :: qmax, enkin, xi, eta, zta, esptar, epsilon 
      logical :: sux 
!-----------------------------------------------
!
!     A routine to split 1 particle in 6
!     Run in 3d kinetic codes
!
 
      fall = fall + 1 
      nsplit = 0 
 
!
!
!     Choose the most energetic particle
!
 
 
      qmax = 0. 
      np = iphead(ijk) 
 6666 continue 
      if (np /= 0) then 
         if (ico(np) == is) then 
            enkin = abs(qpar(np))*(up(np)**2+vp(np)**2+wp(np)**2) 
            if (enkin > qmax) then 
               nsplit = np 
               qmax = enkin 
            endif 
         endif 
         np = link(np) 
         go to 6666 
      endif 
 
 
!
      if (nsplit == 0) then 
         write (66, *) 'splitter failure #1 nsplit=0' 
         return  
      endif 
 
      np = nsplit 
 
      xi = xp(np) - int(xp(np)) 
      eta = yp(np) - int(yp(np)) 
      zta = zp(np) - int(zp(np)) 
      if (int(xp(np))/=i .or. int(yp(np))/=j .or. int(zp(np))/=k) then 
         write (66, *) 'splitter fa' 
         return  
      endif 
      esptar = 1./float(npcel(is))**(1./3.) 
 
      epsilon = min(min(min(min(min(min(xi,eta),zta),esptar),abs(1. - xi)),abs(&
         1. - eta)),abs(1 - zta)) 
      epsilon = epsilon*0.9 
 
      if (epsilon <= esptar/5.) then 
         write (66, *) 'splitter failure: epsilon<0', epsilon 
         return  
      endif 
 
      sux = .TRUE. 
      do ip = 1, 5 
         if (iphead(1) > 0) then 
            npnew(ip) = iphead(1) 
            iphead(1) = link(iphead(1)) 
         else 
            sux = .FALSE. 
         endif 
      end do 
      if (sux) then 
         npnew(6) = np 
!     write(6,*)'splittate',npnew
 
         xpnew(1) = xp(np) + epsilon 
         xpnew(2) = xp(np) 
         xpnew(3) = xp(np) - epsilon 
         xpnew(4) = xp(np) 
         xpnew(5) = xp(np) 
         xpnew(6) = xp(np) 
 
         ypnew(1) = yp(np) 
         ypnew(2) = yp(np) + epsilon 
         ypnew(3) = yp(np) 
         ypnew(4) = yp(np) - epsilon 
         ypnew(5) = yp(np) 
         ypnew(6) = yp(np) 
 
         zpnew(1) = zp(np) 
         zpnew(2) = zp(np) 
         zpnew(3) = zp(np) 
         zpnew(4) = zp(np) 
         zpnew(5) = zp(np) + epsilon 
         zpnew(6) = zp(np) - epsilon 
 
         do ip = 1, 6 
 
            xp(npnew(ip)) = xpnew(ip) 
            yp(npnew(ip)) = ypnew(ip) 
            zp(npnew(ip)) = zpnew(ip) 
            up(npnew(ip)) = up(nsplit) 
            vp(npnew(ip)) = vp(nsplit) 
            wp(npnew(ip)) = wp(nsplit) 
            ico(npnew(ip)) = is 
            qpar(npnew(ip)) = qpar(nsplit)/6. 
!
            xi = xp(npnew(ip)) - int(xp(npnew(ip))) 
            eta = yp(npnew(ip)) - int(yp(npnew(ip))) 
            zta = zp(npnew(ip)) - int(zp(npnew(ip))) 
!
            call weights (xi, eta, zta, wght) 
!
            ipjk = ijk + iwid 
            ipjpk = ijk + iwid + jwid 
            ijpk = ijk + jwid 
!
            ijkp = ijk + kwid 
            ipjkp = ijk + iwid + kwid 
            ijpkp = ijk + jwid + kwid 
            ipjpkp = ijk + iwid + jwid + kwid 
!
            xpf(npnew(ip)) = wght(1)*x(ipjk) + (wght(2)*x(ipjpk)+(wght(3)*x(&
               ijpk)+(wght(4)*x(ijk)+(wght(5)*x(ipjkp)+(wght(6)*x(ipjpkp)+(wght&
               (7)*x(ijpkp)+wght(8)*x(ijkp))))))) 
!
            ypf(npnew(ip)) = wght(1)*y(ipjk) + (wght(2)*y(ipjpk)+(wght(3)*y(&
               ijpk)+(wght(4)*y(ijk)+(wght(5)*y(ipjkp)+(wght(6)*y(ipjpkp)+(wght&
               (7)*y(ijpkp)+wght(8)*y(ijkp))))))) 
!
            zpf(npnew(ip)) = wght(1)*z(ipjk) + (wght(2)*z(ipjpk)+(wght(3)*z(&
               ijpk)+(wght(4)*z(ijk)+(wght(5)*z(ipjkp)+(wght(6)*z(ipjpkp)+(wght&
               (7)*z(ijpkp)+wght(8)*z(ijkp))))))) 
 
         end do 
 
 
         do ip = 1, 5 
 
            link(npnew(ip)) = iphead(ijk) 
            iphead(ijk) = npnew(ip) 
 
         end do 
 
         fall = fall - 1 
         succ = succ + 1 
 
         return  
 
      else 
 
         write (66, *) 'splitter failure: not enough new particles', iphead(1) 
         return  
 
      endif 
      return  
 
      end subroutine splitter 


 
 
 
      subroutine listmkr(np, nploc, link, lvmax, ntmp) 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: np 
      integer , intent(out) :: lvmax 
      integer , intent(in) :: ntmp 
      integer , intent(out) :: nploc(ntmp) 
      integer , intent(in) :: link(0:*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: nps 
!-----------------------------------------------
!
!
!      get particles from linked list &  place indices in table nploc
!
!      INPUTS:
!
!        np      :      first particle index
!        link  :      array of links to next partcile
!        ntmp  :      dimension of nploc array
!
!      OUTPUTS:
!
!        nploc      :      array of particle indices
!        lvmax      :      no of particles in nploc
 
      lvmax = 0 
      nps = np 
      nploc(1) = nps 
    1 continue 
      if (nps > 0) then 
         lvmax = lvmax + 1 
         if (lvmax < ntmp) then 
            nps = link(nps) 
            nploc(lvmax+1) = nps 
            go to 1 
         endif 
      endif 
!
!
      return  
      end subroutine listmkr 
